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ABSTRACT 



We analyze the nonlinear dynamics near the incoherent state in a mean-field 
model of coupled oscillators. The population is described by a Fokker-Planck 
equation for the distribution of phases, and we apply center-manifold reduction 
to obtain the amplitude equations for steady-state and Hopf bifurcation from 
the equilibrium state with a uniform phase distribution. When the population is 
described by a native frequency distribution that is refiection-symmetric about 
zero, the problem has circular symmetry. In the limit of zero extrinsic noise, 
although the critical eigenvalues are embedded in the continuous spectrum, the 
nonlinear coefficients in the amplitude equation remain finite in contrast to the 
singular behavior found in similar instabilities described by the Vlasov-Poisson 
equation. For a bimodal refiection-symmetric distribution, both types of bifurca- 
tion are possible and they coincide at a co dimension-two Takens Bogdanov point. 
The steady-state bifurcation may be supercritical or subcritical and produces a 
time-independent synchronized state. The Hopf bifurcation produces both super- 
critical stable standing waves and supercritical unstable travelling waves. Previ- 
ous work on the Hopf bifurcation in a bimodal population by Bonilla, Neu, and 
Spigler and Okuda and Kuramoto predicted stable travelling waves and stable 
standing waves, respectively. A comparison to these previous calculations shows 
that the prediction of stable travelling waves results from a failure to include all 
unstable modes. 
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1 Introduction 



The dynamics of a collection of weakly coupled limit cycle oscillators can be analyzed 
by treating the oscillator phases {di, . . . , 6n} as an autonomous dynamical system. In the 
context of biological oscillators this approach was advocated by Winfree and subsequently 
developed by Kuramoto for react ion- diffusion systems, [jl], |^ In particular, Kuramoto formu- 



1 



lated the widely studied model for phase dynamics: 

K ^ 

Oi=u;i + -Y.sm{ej-e,)+Ut) (1) 
i=i 

where the natural frequencies {uji} are randomly distributed according to a density g{u!). 
The effect of adding noise C,i(t) to each phase was considered later by Sakaguchi.||] 

For large N, this model describes a transition to collective synchronized behavior as K 
increases above a critical threshold K^. More precisely, in terms of the order parameter, 

1 ^ 

i?(t)e^'^W = -5] e'^^W, (2) 
one finds numerically a transition from an "incoherent" state with R{t) ^ to a partially 



synchronized state with R{t) ~ \/K — K^. for K > Kc. For g{uj) chosen to be unimodal and 
symmetric about u = 0, Kuramoto obtained an analytic expression for the order parameter 
r as a function of the coupling K. His result agrees with numerical simulations and yields a 
simple expression for the critical coupling in the absence of noise 

Kc = (3) 

This transition has been analyzed by several authors and additional references can 
be found in Kuramoto and Nishikawa. If g{uj) has compact support, then for sufficiently 
large K there can be a transition to a completely synchronized state with R ^ l.||12|, |13| 

Since the incoherent state and the partially synchronized state emerge when the initial 
value problem is solved numerically, each must be stable in the appropriate range of K 
at least in an operational sense. However the theoretical explanation of this stability has 
proved to be rather subtle even for the incoherent state, but there has been significant recent 
progress in the work of Strogatz and Mirollo.[^ H Following Sakaguchi, they considered the 
large limit of (|lD and studied the Fokker-Planck equation 

dp djpv) ^ ... 

dt 89 86^ ^ ' 



for the distribution of oscillators. The density p{0, u, t) is defined so that N g{u)p{6, u, t) dO du 
describes the number of oscillators with natural frequencies in [uo^uj + duj] and phases in 
[9, 9 + d9]. Thus p{9, uj, t) d9 denotes the fraction of oscillators with natural frequency uj and 
phase in [9, 9 + d9] and satisfies the normalization 



r d9p{9,u,t) = l. (5) 
Jo 



The identity 



1 ^ 

R{t) sin(7A -9,) = - Y, sm{9, - 9,) (6) 

i=i 

allows the phase velocity of an oscillator to be written as 

v{9,uj,t) = uj + KR{t)sm{iP-9), (7) 
and the order parameter can be re-expressed in terms of the density 

R(^t)e''^m = d9 duop{9,uo,t)g{uj)e''^. (8) 

Jo J~oo 



The diffusion coefficient in (^) reflects the Gaussian noise terms ^i(t) deflned by 

<m> = (9) 

<Us)^M> = 2D6i,6{s-t). (10) 

Equations (^) and (|^ - provide a continuum description of the oscillator population 
for which issues of stability and bifurcation can be analyzed in some detail. The incoherent 
state is described by the uniform distribution po = l/Svr and deflnes an equilibrium for (|) 
since i? = at po- By deflning r], the deviation from po, p{9,uj,t) = po + ^(6*,^^,^) and 
eliminating f in (^), the model can be rewritten as a single equation for rj: 

^ = cv+mv) (11) 



where 



_ d^T] drj Kpo 
Crj = D—r - uj— + 



e-'^ / d9' / diu'T]{9', u', t) g{uo') e'^ + cc 

Jo J-oo 



09^ 09 2 
and the nonhnear terms are 



(12) 
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Note that the normahzation (^) of p imphes 



+e-'' [7] + t^] r d9' r duj'r]{9\ uj\ t) g{uj') e''' \. (13) 



rf^r/(^,w,t) = 0. (14) 



An important quahtative feature of this evolution equation is its symmetry. The group 
0(2) is generated by rotations 

(l)-{9,uj) = {9 + <p,uo) (15) 

and reflections 

K-{9,uj) = -{9,u) (16) 
which act on functions rj{9,uj) in the usual way: 

(7-r7)(^,a;)=77(7-'-(^,^)) 7^0(2). (17) 

Provided g{uj) satisfies 

g{Lo) = g{-u), (18) 

the equation for 77 (|l^) has 0(2) symmetry; more precisely 7 ■ {Crf) = £(7 ■ rj) and '~i ■Mirf) = 
Af{'y-ri) for all 7 G 0(2). A realistic population of oscillators will not have zero mean 
frequency as implied by ([TSD , but if their frequency distribution G(ti;iab) is symmetric about 
the mean ZJ then in the rotating frame defined by 

{9, U) = (6'lab ~Ujt, CJiab - Uj) (19) 



the frequency distribution g{uj) = G{uj + uo) will satisfy (|18|). The densities in the lab frame 
and in the rotating frame are related byQ 

p(^, CU,t) = piab(^ + CJt, CU + Cj,t). (20) 

If g lacks refiectional invariance (|l^) then the symmetry in the rotating frame is reduced to 
the rotations SO(2).[|l3 

Strogatz and Mirollo provide a linear stability analysis of po focussing especially on 
populations with even distributions ([T8| ) that decrease monotonically away from = 0.0 
They show that C has a continuous spectrum with real part equal to —D, and that it may also 
have point spectrum (eigenvalues) depending on the coupling strength. For K sufficiently 
small, the point spectrum is empty, but as the coupling increases a real eigenvalue A emerges 
from the continuous spectrum and for K > K^. moves into the right half plane A > 
signifying linear instability of po. The condition A = determines a critical coupling 



D g{uj) duj 



(21) 



3 D2+^2 

that agrees with results for K^. derived by Sakaguchi using different arguments and reduces to 
(§) when D — s> 0. Thus the onset of synchronization is shown to coincide with the appearance 
of a linearly unstable mode for the Fokker-Planck equilibrium. Subsequently Strogatz et al. 
demonstrated that for D = Q and K < K^, the linearized equation predicts po is stable in 
the sense that R{t) decays via a "phase mixing" mechanism that is qualitatively similar to 
Landau damping in a collisionless plasma. |p 

Recently Bonilla, Neu, and Spigler investigated bifurcations from po assuming a reflection- 
symmetric frequency distribution (^) but relaxing the requirement that g{uj) have a single 
peak.|TI|] For bimodal distributions they showed that C could have complex eigenvalues so 



^In the lab frame the evolution equation for piab(^iab, t^iab,i) has rotation symmetry S0(2) but not 
reflection symmetry. Calculations done using lab frame variables have puzzling features that are non-generic 
for SO (2) symmetry and have their origin in the full 0(2) symmetry of the rotating frame. 



that the instabihty of po occurs via Hopf bifurcation. Amphtude equations were derived to 
describe the nonhnear development of the instabihties, but unfortunately the Hopf analysis 
retains only half of the unstable eigenvectors and as a result nonlinear travelling and stand- 
ing wave solutions were overlooked. The fact that the 0(2) symmetry leads to degenerate 
eigenvalues of multiplicity-two was apparently not recognized. 

Independently, Okuda and Kuramoto analyzed the onset of collective behavior in a pop- 
ulation comprised of two sub-populations using a different model; each sub-population had a 
unique frequency in the absence of coupling and was described by a number density. The 
instability of the incoherent state was studied under several assumptions on the couplings 
between the oscillators. In the regime where the coupling between the populations equals 
the coupling within a population, their model is a special case of (^). They also find that 
instability can arise from either real or complex eigenvalues and derived amplitude equa- 
tions for the nonlinear behavior. Their equations include both travelling and standing wave 
solutions although the role of symmetry in the problem is not discussed. 

In this paper I re-examine the bifurcation problem and derive amplitude equations for 
both steady state and Hopf bifurcation using center manifold reduction. In Section II, the 
spectrum of L is discussed to clarify the treatment of complex eigenvalues and to discuss the 
occurrence of non-semi-simple real eigenvalues. These arise naturally at the multi-critical 
point where the Hopf and steady-state bifurcations coalesce. This codimension-two point 
corresponds to a so-called Takens-Bogdanov bifurcation and the relevant amplitude equa- 
tions will be investigated in future work. The discussion of Section II is facilitated by the 
similarities of L to the linear operator for the Vlasov-Poisson equation which has been dis- 



cussed elsewhere. |15 



In Section HI, the amplitude equations are derived and their limiting form as i5 — 
determined. The results are strikingly different from analogous amplitude equations for 
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Vlasov instabilities; in the collisionless plasma case the coefficients diverge as the critical 
eigenvalues approach the continuum but in this model the coefficients have finite limits. ||T6[ 
In the nonlinear theory, the presence of 0(2) symmetry is exploited to take advantage of 
the well-established theory of bifurcation in systems with circular symmetry. |T^, |TB[ For a 



refiection-symmetric bimodal population comprised of two Lorentzians, we evaluate the cubic 
normal form coefficients. For the steady-state bifurcation the appearance of the synchronized 
state can be either supercritical or subcritical. In the Hopf bifurcation we find that the 
standing waves are supercritical and stable while the travelling waves are supercritical and 
unstable. These conclusions are compared to previous studies. [|], |ll 



2 Linear theory 

The Fourier expansion of rj, 

oo 

V{e,uj)= Y: Vi{uj)e''' (22) 

(=— oo 

allows the linear operator ([T2|) to be re-expressed as 

oo 

E i^iVi)e'^' (23) 

i=— oo 

where 

/oo 
duj'g{uj')rii{uj'). (24) 
-oo 

Note that the normalization condition (|T^ on rj implies 

Vi=o = 0. (25) 
An adjoint operator (A, CB) = {C^A^ B) can be defined for the inner product 

p2'K roo 

{A,B)= dO du A{e,uoy B{e,u). (26) 



In terms of the Fourier expansion A = J2 ^^e*'^, this definition yields 

oo 

C^A= (LM/)e*'' (27) 



l=—oo 

where 



oo 



(lJA)H =a{u + ilD) Ai{uj) + ^^^{6i,i + r~ du'Aiiuj'). {2t 



It is convenient to introduce a concise notation for the integration over u in pq ) 

roo 

<A,B>= / dujA*B; (29) 



oo 



then and L| satisfy < A, LiB >=< LjA, B > . 
The eigenvalue problem for C, 

/:^ = A^, (30) 
can be treated on each Fourier subspace separately so we set 

^(^,cj) = e^'VM, (31) 

and reduce (|30|) to 



Liij = Xip. (32) 
After scaling the eigenvalue by A = —il z, the eigenvalue equation (^) becomes 

[{u -z)- ilD] ^ = ^— (5,,! + (5i,_i) / du'giu') ^{u'). (33) 



The first observation is that |/| = 1 is a necessary condition for to be a well-behaved 
eigenfunction. When |/| ^ 1 then (33) becomes [{uj — z) — ilD] ip = 0, and we have ip = 



unless Im^; = —ID. For Imz = —ID then [{uj — z) — ilD] is non-zero except ai u = Kez, so 
the only solutions are delta functions 6 {u — Re z) . These solutions correspond to a continuous 
spectrum for along ReA = / Imz = —l^D; as D ^ the continuous spectrum moves onto 
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the imaginary axis. For |/| 7^ 1 this is the only contribution to the spectrum of L and it does 
not depend on the couphng K or on the frequency distribution g{uj). 

For |/| = 1 it is possible to have eigenf unctions that are non-singular; their occurrence 
does depend on K and g{u]). Now the eigenvalue equation (|55|) becomes 



— 7/ K 

[{uj-z)-UD]^ = ^ du'g{J)^{u'). (|/| = 1) (34) 

Z J —00 

When Im^ = —ID, then as before we get singular solutions corresponding to a continuous 
spectrum along ReA = —D\ the detailed form of these solutions will not be needed. If 
Imz ^ —ID then {uj — z) — ilD cannot vanish, and we have 

If the integral vanishes then ip = Q so we may assume the integral is non-zero and normalize 
iIj such that 

/oo 
du'g{u')ip{u') = 1. (36) 
-00 

Consistency between (BB) and (BBI) requires 2; to be a root of the "spectral function" 



A,(,)^l^^lEr ,U^^. = (37) 

2 J -00 UJ — z — tlD 

The spectral function A; satisfies two identities; for arbitrary g{u) we have 

A^izY =A_iiz*) (38) 

so when (2:0) = determines an eigenvalue A = —Uzq then A_i (zq) = determines the 
complex-conjugate eigenvalue A*. In addition when the frequency distribution is reflection- 
symmetric (0), the problem has 0(2) symmetry and A; satisfies 

Aiizy = Aii-z*). (39) 
Now roots come in pairs whenever Re z 7^ 0, i.e. whenever A is complex. 



The roots of (z) determine the eigenvalues and the nature of the root determines 
whether the eigenvalue is semisimple or not. Recall that the projection operator Pa onto 
the generalized eigenspace Ex for an isolated eigenvalue A = —Uzq is defined by a contour 
integral,]!^ 

PxA = ^ f d\'{X' -C)-'A- (40) 

where Fq is a small loop enclosing A in a counter-clockwise sense. The resolvent operator 
(A — £)^^ can be expressed a^ 

oo 

iX-C)-'A= Y: (RiiX) Ai) e^^' (41) 

Z=— oo 

where -R;(A) = (A — L^)"^ denotes the resolvent for Li and can be calculated explicitly 
n(X)A I i\i + 5i,^i)K/2 .00 giuj')Aiuj') 



With ( PI) P\ may be simplified to 

id —id 

-hi ''^'«'(^''-^' ^hL (43) 

since eigenvalues can only occur for |/| = 1. If -Ea, the range of Pa, is one-dimensional then 
A is simple; when the dimension is greater than one and the eigenvectors of C span E\ the 
eigenvalue is semisimple. If dimE'A > 1, and the eigenvectors do not span Ex, then A is 
non-semisimple. Simple roots of A; {z) lead to simple or semisimple eigenvalues and multiple 
roots give eigenvalues that are non-semisimple. We discuss this point in detail for the case 
of a reflection-symmetric frequency distribution and 0(2) symmetry. 

2.1 Simple roots and semisimple eigenvalues 

For a simple root zq of Ai {z), we have 

Ai (zq) = a; (^0) ^ (44) 

^In this sum, the term I — may be omitted since the relevant space of fmictions satisfies (p 

10 



where A' denotes the derivative of Ai (zq) with respect to z. The identities (pSD - ( |39D imply 
the additional simple roots 

Ai i-z^) = A_i (z;) = A_i (-^o) = 0. (45) 

If Re 2:0 = (a real eigenvalue) then zq and —Zq are the same root and roots {zq, —zq) also 
coincide. When Re 20 7^ (complex eigenvalues) then these roots are all distinct. 

2.1.1 Real eigenvalues (Re = 0) 

With Re Zq = 0, the roots zq = and Zq for / = ±1 lead to linearly independent eigenvectors 
for the real eigenvalue A = —izQ. From (^) and (^^, these vectors are 

^{e,uj) = e'^^iuj), (46) 

where 

^(^) = , -n (47) 

[UJ — Zq) — tJJ 

and the conjugate function 

*(^,cj)* = e-*V(^)*- (48) 

Note, that these eigenvectors are related by reflection '^{6, uj)* = {k-'^){9, u). Corresponding 
to \1' and K, ■ \E', there are adjoint eigenvectors 

C^^ = (49) 
C\k-^) = A(k-^) (50) 



where 



with 



^{9,uj) = ^^P{u;) (51) 
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These eigenvectors are conjugate k ■ ^ = ^* and satisfy 



K ■ \1>) 



(k ■ ^, K • ^) = 1 
(k- ^,^) = 0. 



(53) 

(54) 



For example, the eigenvalue equation ( ^91 ) can be verified directly from the definitions ( pT]) 

- m 



2vrA'i (2;o)* 



ig{uj) 



{uj-z^ + iD) 



g{u)K 



271 k'^ {z^y 



2 J —CO 

-ig{uj) 
ig{uj)zl 



du'- 



uj' - + ^D) 



{uj-z^ + iD) 



tgiu)iA.^iz*)-l) 



iz, 



-9(^)4 



° 27rA'i(^o)* \{uJ-z^ + iD) 



(55) 



since A = A* = —izQ. 

Since Ai (2:0) = corresponds to a simple pole, the result of projecting rj onto is easily 
evaluated. From (|42|) - (|^) we have 



iO 10 



Ad 



dX' 



K/2 



2m Jfo Ai {iX') {icu + D + A') J-oo iu' + D + X' 



2Tii Jfo A_i{-iX'){-iuj + D + X') J-oo -ito' + D + X' 
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K/2 



iA'i {zo){iuj + D — izo) J-oo iuo' + D — izo 

-zA'_i {z^){-iiu + D + iz^) J-oo -iuj' + D + iz*^ 



= {^,r])^{6,u) + {K-^,r]) K-'^{9,uj). (56) 

The eigenvectors \1/ and k ■ ^' clearly span the range of Pa- 

2.1.2 Complex eigenvalues (Re 20 7^ 0) 

When Re^o 7^ 0, then zq and —Zq are distinct roots of Ai (z) and A_i (z), respectively, 
that determine two linearly independent eigenvectors for the complex eigenvalue A = —izQ. 
These eigenvectors are \E' given by (^61) - ( ^Tf ) and the reflected function 

K-"^ = e-'^-^i-uj). (57) 



This second eigenvector was omitted in the Hopf analysis of Bonilla et al. |Tl[| The remaining 
roots —Zq and Zq for Ai and A_i determine the conjugate eigenvectors ^* and (k ■ \I')* 
corresponding to the conjugate eigenvalue A* = iz^. The adjoint eigenvectors \E' and k ■ ^ 



have the same form (^) - ([54 ), however now (|49|) - (|50D become 

£t(^.^) = A*(k-^); 



(58) 
(59) 



the adjoint eigenvectors for A are the conjugate functions and (/«■ \E')*. The orthogonality 
i3) remain valid. 



relations (|53D - ( 

The projection PxTj can be evaluated as before (^) except that now the relevant root 
(pole) for Z = — 1 is A_i {—zq) = 0. The final result has the same form 



P^T] = {"^ , 1]) + {k ■ , 1]) K ■ 



(60) 
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A similar evaluation of Pa* V gives 

Px*V = {PxVy- (61) 
In each case the eigenvalue is semisimple with multiplicity- two. 
2.2 Multiple roots and non-semisimple eigenvalues 

A single root, Ai (zq) = 0, determines exactly one eigenvector p6| ) - ( |T7| ) for the eigenvalue 
A = —izo- There can be additional linearly independent eigenvectors for the same eigenvalue 
A only if A_i simultaneously has a root at —zq. In any event, a given eigenvalue can never 
have more than two eigenvectors. However if the root zq of Ai is not simple, the range of 
Pa turns out to have dimension greater than the number of eigenvectors because there are 
generalized eigenvectors in addition to the true eigenvectors. A generalized eigenvector T 
satisfies 

{C-Xyr = (62) 
{C-Xy-^T (63) 

for some integer j > 2. Note that ( |62D simply says that {C — A)''^^T is an eigenvector for A. 

Since degenerate roots of Ai arise when simple roots collide, we expect generalized eigen- 
vectors to enter the linear problem whenever eigenvalues collide as parameters are varied. 
The connection between degenerate roots and generalized eigenvalues can be made explicit 
as follows. Consider the simplest situation j = 2 and let \E' = e^^ip be a true eigenvector 
corresponding to Ai (^o) = 0. We seek a generalized eigenvector T = e^^v satisfying 

{C-X)T = ^ (64) 
{C-XfT = 0; (65) 
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note that ( |6^ implies (]65|). From the definition of £, (|6^ ) can be rewritten as 

{u-z^- iD)v + — / duj'g{uo') v{uj') = (66) 



upon multiplying by g{uj)/{uj — Zq — iD) and integrating ( |5BD becomes 

/■oo 

Ai (^o) / duj'g{uj') v{uj') = -tA[ (zo). (67) 



Since Ai (zq) = 0, the assumption (|6^ ) that a generalized eigenvector T exists leads imme- 
diately to a contradiction if the root is simple, i.e. if A'^ {zq) ^ 0. An extension of this 
argument shows that a generalized eigenvector with j >2 for A cannot occur unless the first 
j — 1 derivatives of Ai vanish at the root. 

Bonilla et al. found that the transition between a bimodal distribution g{uj) supporting 
real eigenvalues and a bimodal distribution g{uj) supporting complex eigenvalues occurs via 
a non-semisimple real eigenvalue corresponding to a double root: 

Ai(^o) = A;(zo) = Ai (^0)7^0 (68) 

where Re 2:0 = 0. For more complicated distributions, it is likely that double roots with 
Re 2:0 7^ are also possible. The case of a real eigenvalue is considered here and we assume 
Re 2:9 = for the remainder of this section. The double root (|68D at zq then implies a double 
root at —Zq for A_i: 

A-i(-^o)=0 A'„i(-^o) = A':i(-zo)^0. (69) 

In applications with double roots, the frequency distribution will typically depend on at 
least one additional parameter cuq, for example in a bimodal distribution uoq could be the 
distance between peaks in g{ijj^ujQ). Then A'^ = determines zq as a function of K and 
ujq and Ai = implies a relation K = K{ujo) between the parameters. The eigenvector 
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corresponding to Ai (^o) = is then \1/ = e ip with 



(70) 



{u! — Zq) — iD ' 

and the generahzed eigenvector T = e^^v is determined by 

^(^) = 7 ('^1) 

(a; — 2;o — xL'j^ 

One can always add to v{uj) an arbitrary multiple of the homogeneous solution ip but without 
loss of generality we have set the coefficient of this term to zero. Reflection by k, gives the 
eigenvector and generalized eigenvector corresponding to the second root (^): k ■ = 
e~^^ ■ip{—iu) and k-T = e~^^ v{—uj). 

For the adjoint eigenvalue problem, there is an eigenvector (£^ — A) \& = and generalized 
eigenvector {C^^ — A) T = \E' given by \E' = e"^^ ip{uj) and T = e^^ v{uj) where 

-2g{uj,uJo) , 2A'l'{zo)g{uj,uJo) 



A'l {zo){uj - ^0 - toy ^ 3A'/ {zoYiu ~ z^ - iD) 
2ig{uj,uJo) 



(72) 
(73) 



A'({zo)iiu-Zo-tD)_ 

In addition the reflected vectors k ■ \E' and k, ■ T solve the same equations. These eigenvectors 
satisfy the biorthogonality relations 



(^, V[/) =< ^ >= 1 =< K- ljj,K-tp >= 1 

(T ,T) =< v,v >= 1 (k-T, K-T) =< K ■ V, K ■ V >= 1 



(74) 
(75) 



and 



{^,T) =<ip,v >=0 (k-^,k-T) = 
(f,^) =<d,ip >=0 (/t-t,K-*) = 0. 



(76) 
(77) 
(78) 



16 



The evaluation of Pa V is now more involved because the poles are not simple; from 
we have 



2TTi Jro Ai {iX') {iu + D + A') J^oo iu' + D + X' 

+ I dX' r dj 9{^'^^^)^-^(^') (rg) 

2m Jro A-i {-iX') {-iuj + D + X') J-oo -iuj' + D + X' ' ^ ' 

Consider the first integral for K = K{ujo) when Ai has a double root; the contribution from 
this double pole is 



2m Jro Ai {iX') {iu + D + X') J-oo iu' + D + X' ^ ^ 



dX' 



iX'-X)^KM/2 \ r ^^,9iu;\uJo)vi{^') 



Ai{iX'){iuj + D + X') [J-oo iuj' + D + X' 



^1) 



J A'=A 



Expanding Ai, 



Ai (^A') = t^A'l {tX) ^ + ^^^A;" (^A) ^ + ■ ■ • , (82) 



3! 

and evaluating the right-hand side of (|8TD gives 



2m J To Ai[iX'){iuj + D + X') J-oo lu' + D + X' 

A similar evaluation of the second integral in (^) yields {^,7])'^ + {T,7])T which is the 
complex conjugate of (E^. Overall we have 



Pa ^ = r/)* + (T, r/)T + (^, r/)^ + (T, r7)T; (84) 

thus a purely imaginary double root determines a real eigenvalue with two eigenvectors and 
two generalized eigenvectors. 
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Complete expansions in terms of the eigenfunctions of L/ and hence L can be derived as in 
the Vlasov case[|l5| if desired. We do not need these expansions for the bifurcation analysis, 
it is sufficient to know how to project out the components of r/ along the critical eigenvectors 
in the three cases of interest: semisimple real, semisimple complex, and non-semisimple real. 

2.3 Example of a bimodal population 



We illustrate this discussion for the population 

1 



5((u;,u;o) = ^ 



+ 



1 



5) 



comprised of two sub-populations centered at = icjo; the e — limit yields the delta 
function distribution considered independently by Okuda and KuramotoQ and by Bonilla 
ei a/.|TT[ When cJq > there are two peaks in giuj^uj^^ located at 



2- + 



(86) 



For Im z > —D, the evaluation of Ai yields^ 

[z + z(e + D)f - _ ^ ^ ^)]/2 



Ai W 



which has two roots for ujq > 0: 



z± 



[z + i{e + D)f-ujl 



(87) 



[K/A)-e-D±J{K/AY-ujl 



At ujq = 0, the expression in ( |87[ ) has only a single root z+ = i [{K/2) — e — D]; the solution 
z_ at —i (e + D) is an artifact. For < 4ujq < K, the roots are pure imaginary and the 
eigenvalues are real; along the line 4ujq = K the roots collide, a double root forms, and a 
non-semisimple real eigenvalue occurs. When 4ujq > K, the eigenvalues are complex and 

Hopf bifurcation is possible. 

■^This case is easily treated by closing the integration contour in the lower half plane. One can show that 
for Im z < —D and K > 0, the spectral function has no roots. 



For Aujq < K, the root in (^) determines an eigenvalue if Im Zj^ > ~D which 
requires either K > 2 {uJq + e"^) / e or K > 4e; the solution z^ determines an eigenvalue when 
ie < K < 2 (ciJq + e^)/e. The eigenvalues X± = —i z± disappear into the continuous spectrum 
at — D when A± = —D. This defines a boundary in parameter space 

K = Ke = 2{ujl + e^)/e (89) 

which is tangent to = Aujq &i uoq = e + D . The real eigenvalues cross the imaginary axis 
A± = along a second boundary 

K = K, = 2 [ujI + (e^ + D^)]/{e + D) (90) 

with A_|_ = if < 4(e + D) and A_ = if > 4(e + D). This second boundary is tangent 
to K = Aujq dX K = 4(e + D) where we find a non-semisimple eigenvalue at zero; here the 
steady state bifurcation turns into Hopf bifurcation. Figure 1 shows the {K, ujq) parameter 
space with the curves corresponding to eigenvalues coalescing with the continuum, crossing 
the imaginary axis, or colliding on the real axis. 



For Aujq > K, the roots ([88| ) are complex and the eigenvalues (A, A*) are given by 

A = -izo = {K- Kh)/A + i n (91) 



where Q = — (K/A)"^ and Kh = 4(e + -D). If K decreases below K'^ = 4e the eigenvalues 
(A, A*) disappear into the continuous spectrum at Re A = —D; as K increases above Kh 
there is a Hopf bifurcation. The curves K = K'^ and K = Kh are also shown in Fig. 1. 

In an early discussion of bimodal populations without extrinsic noise D = 0, Kuramoto 
suggested that once the maximum of g{uj) was no longer at the average frequency (here 
uj = 0) the incoherent state would go unstable for some K = K^q < K^ with synchronized 
behavior nucleating among oscillators whose frequencies lay near the maxima of (7(^7). 
Although the development of bimodal peaks in ( ^5]) for uq > y/e/3, does alter the stability 
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properties of the incoherent state; the effects are somewhat more comphcated. As we shall 
see below in our analysis of the steady state bifurcation, for 



the incoherent state is stable when K < and there is a subcritical steady state bifurcation 
to a synchronized state at = Kc- There is no prior onset of synchronization in the 
peaks of the distribution. The further qualitative change in the stability properties of the 
incoherent state with increasing uq above (e + D) is to replace the steady state bifurcation 
by a Hopf bifurcation which leads to time dependent behavior (cf. Section 3.2) in the form 
of stable nonlinear "standing waves" . These waves appear to be the type of state Kuramoto 
envisioned; they are discussed in greater detail below. It seems likely that these conclusions 
remain true in the limit D ^ 0. 

3 Codimension-one bifurcations 

For the nonlinear analysis to follow we require the Fourier expansion of A/'(?7). With 




(92) 




(93) 



and 




— oo 



oo 



duj'ri{9', Lj', t) g{uj') e 



2n < g,r]i >, 



(94) 



Af{ri) in (|T3|) becomes 




+ ^ e*'^ [< 9. V^i > Vi+ii^)- <g,Vi> Vi-ii^)]}- 

\l\>2 



(95) 
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3.1 Steady-state bifurcation 

We consider a simple pure imaginary root Ai (2:0) = so that A = —izo is real. As 
a parameter such as K is varied, we assume that A crosses through X = at K = K^. 



From (^q ) - (^SQ there are two eigenfunctions ^ and k ■ \E' = \1'*, with amplitudes (a, a*) 
defined by a{t) = rj). The distribution 77 decomposes into the critical linear modes Px rj = 
a{t)'^{9, uj) + cc, and the remaining degrees of freedom S{9, u) = r] — Px rj. Thus we have 

ri{e, u, t) = [a{t)<iJ{9, uj) + cc] + S{9, u, t) (96) 

where (\^, S) = 0. In terms of the Fourier components of rj and S, we have 

r]i = ai!6i^i + a* ip* 6i^^i + Si; (97) 

note that the normalization (|25|) requires 

^/=o = 0. (98) 

By projecting with the evolution equation (^) can be rewritten as 

a = Xa+ {■^,^^{r|)) (99) 
— = £S + {^/i7^)~[{^,A^{v))^ + cc]}; (100) 

since < g,rii >= a 61^1 + a* < g,Si > in (p5|), the Fourier decomposition of these 

equations gives 



a = Xa-TcK {a*+ < g,S^i>) <iIj,S2> (101) 
S 
dt 



(99 

— = CS-7rK{e'\a*+<g,S-i>){S2-<i',S2>^) 



+2 [{a*+ < g, S_i >)53(^) - (a+ < g, 5i >)(« ^ + 5i)] 
+ ^ /e*'^ [(«*+ < >)5;+i(^) - (a+ < >)Si^i{u;)] 

l>3 

+cc} . (102) 



21 



3.1.1 Dynamics on the center manifold 



The equations (|101|) -( p!02l) are obviously coupled but on the center manifold the "trans- 



verse" degrees of freedom S depend on time through the critical amplitudes; this may be 
expressed in terms of a function H{9, u, a, a*): 

oo 

S{9, u, t) = H{e, to, a{t), a{ty) = ^ Hi{9, a{t), e*'^ (103) 

In geometric terms, H describes the graph of the center manifold near the incoherent 
state. ||20|, ^ Thus solutions on the center manifold satisfy (|1Q1| ) - (|102| ) with S replaced 



by if 

a = Xa-'KK{a*+ < g,H_i>) <ijj,H2> (104) 
— = CH-7rK{e'%a*+<g,H^i>){H2-<i^,H2>tlj) 

+2 e'^' [{a*+ < g, >)H,{u) - {a+ < g, H, >){a^ + H,)] 
+ J2 I [(«*+ < 9. ^-1 >)^/+i(^) - («+ < 9. Hi >)Hi^i{uj)] 

l>3 

+cc} ; (105) 
the equation ( |104| ) for the mode amplitude now decouples from 5* and defines an autonomous 



two-dimensional system. 

The center manifold dynamics (|104| ) and the graph function are both constrained by the 
0(2) symmetry of the problem. The Fourier coefficients of H for / > must have the form 

Hi{uj, a, a*) = a'' hi{uj, a); (106) 

where hi{uj,a) satisfies 

hi{-u;,ay = hiiu;,a) (107) 

and a = jap denotes the basic 0(2) invariant. Similarly ( |104|) , the equation for a must have 
the general form 

a = p{a)a (108) 
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where p{cr) is a real-valued function of a single variable. For small a, the leading terms in 



|^ |104| ) are therefore 

a = a[p{0) + p^iO)\a\^ + ■ ■ ■] (109) 



and comparison to (|104|) shows p{0) = X and Po-(O) remains to be calculated. 

We briefly summarize the derivation of (|106|) - (|108|) . The form of H is constrained by the 
fact that the center manifold is mapped to itself 0(2) transformations, [p^j More precisely, if 
ri{6,uj) corresponds to a point on the center manifold and 7 G 0(2), then the transformed 
function 7 ■ rj will also lie on the center manifold. When we represent elements of the center 
manifold using ( p.03| ). 



ri{9, iu) = a ^(^, to) + a* ^*(^, cu) + H{e, to, a, a*), (110) 

then 7 ■ 7] must also be of this form 

(7 ■ 7]){e, 00) = a' ^{6, 00) + a'* ^*{e, uj) + H{e, 00, a', a'*) (111) 

for some transformed amplitude a'. The constraints on H implied by this requirement follow 
from considering the generators of 0(2) in (^) - (16). In conjunction with ( |T7|) we find for 
rotations 

[a ^{6 - 0, cj) + cc] + H{e - 0, 00, a, a) = [a e''^ ^(^, u) + cc] + H{e, e''"^ a, e'^ a) (112) 
and 

[a *(-0, -to) + cc] + H{-9, -u, a, a*) = [a* ^(^, to) + cc] + H{9, to, a*, a) (113) 

for reflection. Since \E'(6' - 0,cj) = e~''^'^{9,uj) and "^{-9, -u) = '^{9,u)*, the amplitude a 
transforms under 0(2) by 

(f).a = e-'^a (114) 
n ■ a = a*, (115) 
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and H must satisfy 

H{e -(f),uj,a,a*) = H{e,uj,e-''^a,e'^a*) (116) 
H{-e,-cu,a,a*) = H{e,uj, a* , a). (117) 
Restated for the Fourier coefficients Hi, these relations require 

e-'^^Hi{u,a,a*) = Hi{e,uj,e-"l' a,e'^ a*) (118) 

and 

Hi{-uj,a,a*y = Hi{uj,a*,a), (119) 

respectively. The rotational constraint ( |118| ) states that the function {a*y Hi is invariant 
under rotations, hence it depends on the mode amplitudes {a, a*) only through the basic 
invariant^ cr = |ap and must contain an overall factor of (a*)': 

{a*yHi{uj,a,a*)=a^hi{uj,a) (120) 

where hi{uj, a) is an arbitrary function of u and cr. Dividing by («*)^ gives the form of Hi in 
( |106|) ; the reflection symmetry (|119|) then further requires that hi{uj,cr) satisfy (|107|) . 



The 0(2) invariance of the center manifold also implies that the vector field in ( |104| ) 
must commute with the representation of 0(2) given in ( |114| ) - ( |115| ). This can be explicitly 
verified using the properties of Hi in (|118|) - (|119|) . Any smooth two-dimensional vector 



field that commutes with this representation must have the form ( |108| ); see for example 
Golubitsky, Stewart, and Schaeffer, chapter XII, §5. 

3.1.2 Evaluation of p^(0) 



From (|106|) , Hi and H2 have the form 

Hi = ahi{uj,a) =a[hf\uj) + 0{a)] (121) 
H2 = a' hiuj, a) =a^h^^\u)+0{a)], (122) 



*cf. Golubitsky, Stewart, and Schaeffer, chapter XII, §4 
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however the geometry of the center manifold requires that the lowest order terms be quadratic 
in {a, a*) so hf'^ = 0; hence the center manifold vector field ( |104| ) becomes 

a = Xa-'KKa\a\^{l+ < g,hl>) <ijj,h2> (123) 
= a[X-7TK <i),hf > lal"^ + ■■■]. (124) 



Comparing (|109|) and ( |124| ) we obtain Pa{0) 

pM = -ttK <^,h^? > (125) 

(2) 

in terms of the center manifold coefficient ■ 
When S is given by ( |103| ) we have 

dS dH dH 



for solutions on the center manifold, and consistency with (|102|) requires 

—a + ^a* = CH-7cK\e'%a*+<g,H_i>){H2-<i!,H2>^) (127) 
oa oa* ^ 

+2 e*2e < H^, >)H,{u) - {a+ < g, H, >) (a ^ + H^)] 

l>3 

+cc} 



where a is given by (|104|) . This provides the equation needed to determine the expansion 



coefficients of H; for small |a|, we can solve (|127|) using a power series expansion in (a, a*). 
Separating the different Fourier components of ( |12?| ) for 1 = 1,1 = 2 and / > 3, we find 

FiTJ PlTJ 

^d + ^d* = UHi-TxK{a*+ <g,H_i>){H2- <^,H2>i^) (128) 



^d + ^d* = UH2-2TiK[{a*+<g,H_,>)H^{u) 
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-{a+ <g,Hi>){ai) + Hi)\ (129) 



-ia+ <g,Hi>)Hi^i{uj)] (130) 



respectively, and inspecting these equations shows that the inhomogeneous term a'^ip in (|129|) 
that forces a non-zero term h^2\^) ^he lowest order term in H2. Since hf\uj) = 0, the 
balance at second order in (|129|) is 2A h'^^ = L2h^2^ + 271 K ip which gives 



(2), . -mKilj{uj) 
^2 M = 131 



Thus the cubic coefficient in the center manifold dynamics (|125|) is given by 



= , , / duo- -. 132 

' 2 A'l (zo) i-00 {uo-zo- iDf {uo-zq- i2D) ^ ' 

Note that Po-(O) is real as expected from 0(2) symmetry. 
3.1.3 Example of a bimodal population 



For the population in ( poD the first steady state bifurcation occurs when the eigenvalue 
A+crosses zero; this occurs along K = Kc Fig. 1. Along this steady state locus, the coefficient 
( |132|) can be easily evaluated by closing the contour in the lower half plane: 

2[(e + DY-uMe + 2DY + u:l] ^^''^ 

where 



Since c^o < (e + D) when this bifurcation occurs, the overall sign of Po-(O) is determined by 
uJq — uj'^. For ujq < ujdi the bifurcation is supercritical to a stable synchronized state: for 
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ujo > ujd, we have Pcr{0) > and the bifurcation is subcriticaL Note that Pa{0) is finite in 
the hmit D ^ as long as e > 0; this feature is examined in greater generahty below. 

For e > and D = 0, the transition from supercritical to subcritical bifurcation was 
noted by Kuramoto,0 and at e = and D > it was found by Bonilla et al. and Okuda 
and KuramotoJlTTI, The former group also calculated the fifth order term necessary to 
determine that the subcritical branch turns around at a saddle node bifurcation and becomes 
stable at finite amplitude. 

3.1.4 The D ^ limit 

If D = then the critical eigenvalues emerge from the continuous spectrum simultane- 
ously with the onset of the linear instability, and it is of interest to examine the behaviour 
of the amplitude equation in this regime. At criticality for either steady-state or Hopf bi- 
furcation ReA = so Imzo = 0, thus we consider the behavior of P(t(0) when D Q with 
Zq = 0. The integrand in (|132| ) has poles at uj = iD and uj = i2D that approach the contour 



of integration as D ^ 0. We assume g{uj) is analytic at = so that the contour can be 
deformed slightly below the axis and the limits evaluated using the Plemej formulas. 

It is convenient to first consider the limiting behavior of A^. From (^) we re-express the 
derivatives of A^ (z) by integrating by parts 

A; (z) = 5n,o + — / duj ^ • (135) 



(136) 



dz^ ' 2 {uj-z-iWy 

thus when 2; = r is real-valued, then as D ^ 



lim - — Al (r) = Ono H 



00 ijj — r 010"- 



where V denotes a principal value integral. For steady-state bifurcation, when g{u) is even 
in uj and r = at = K^., these limits simplify: 

„™ |1a, (0) = r „ odd (137) 

D^Odz'' ^ ^ Jo UJ ^ ^ 
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lim (0) = 6no + n even. (138) 

Note that for = by setting Ai (0) = in ( |138| ) recovers the previous result (^) for 
K, = 2lT,g{Q). 

For small D the behavior of Po-(O) is extracted in the same way; setting 2:0 = we find 



^ 2 A i3-0i_oo 



2A'i(0) D~.oJ^oo {00 - iOy {io - i2D) 



.2 A'/ (0) 



z (TTiT,)" + o(D) 

^ ' 2Ai(0) ^ ' 



-TT 92^7(0) 



^(0)2 



°° , dgiuj) duj 
doj 

U) 



+ 0{D). (139) 



Once differences in notation are taken into account, this result agrees with the coefficient 
obtained by Kuramoto by the "self-consistent" method. |@] Notably absent is the kind of 
singular behavior found in the Vlasov amplitude equations in the corresponding limit. |]16| For 



a monotonic reflection-symmetric profile this gives P(j(0) < in the limit of weak diffusion. 
Thus at the level of the amplitude equation, we expect supercritical bifurcation to a stable 
synchronized state in this regime. However when D = Q center manifold theory no longer 
justifies our reduction to two dimensions; the qualitative agreement ai D = Q between 
numerical simulations 0] and our amplitude equation may be fortuitous. 

3.1.5 Perturbing 0(2) S0(2) 

If the distribution g{oj) is perturbed so that the reflection symmetry (|T^) is broken, then 
the bifurcation problem (in the rotating frame (^)) has only the rotational symmetry SO (2). 
In this case one expects the real eigenvalue of multiplicity-two to split into a non-degenerate 
complex conjugate pair and the steady state bifurcation considered here to be perturbed 
to a Hopf bifurcation. The resulting Hopf bifurcation leads to time-periodic states in the 
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form of rotating or travelling waves; Knobloch discusses such perturbed bifurcations in more 
detail. 



3.2 Hopf bifurcation 

We next consider a simple complex root Ai (zq) = so that zq ^ —z*^ and A = —iz^ 
is complex. As K is varied through K^^we assume that Re A crosses through ReA = 
aX K = Kh- From (^) there are two eigenfunctions \1' and /t ■ \E' = \E'*, with amplitudes 
defined by a{t) = {^,r]) and f3{t) = {k ■ ^,r]). With the projections (H) - (|TD, the 
distribution r] decomposes as before, 



riie, cu, t) = [a{t) ^(^, iu) + p{t) K ■ ^{9, cu) + cc] + 8(6, cu, t) 



where S) = {k ■ S) = 0. 

Projecting ([TI|) with ^ and k ■ ^ yields 



a 



dS 



In terms of the Fourier components of rji 

r]i = {a^ + (3* K-^*)Si,i + {a*ilj* + (3k-iIj)Si,.i + Si 
<g,m> = ia + /3*)6i,i + {a* + f3)6i,-i+ < g,Si>, 

and the expansion oiAfirf) in (^), ( |141| )- (|143|) can be re- written as 



a 

dS_ 
'dt 



(140) 



(141) 
(142) 
(143) 



(144) 
(145) 



Xa - ttK {a* + [3+ < g, S^i >) < ^, 5*2 > 
Xp -7rK{a + (3*+ < g, Si >) < k ■ tp, S^2 > 

CS - TvK [e'\a* + (3+ < g,S_i >)[S2- < ^^,^2 > tp- <k-^%S2 > K-tfj*] 



(146) 
(147) 



29 



+2e'^'[{a* + 13+ < g,S-i >)S3H - (« + /?*+ < g,Si >){a^ + P* k ■ ^ + Si)] 

l>3 

+cc}. (148) 

3.2.1 Dynamics on the center manifold 

For a solution on the center manifold near the incoherent state, we can express the 
time-dependence of S{9,u!,t) as 

oo 

S{e,u;,t) = H{9,uj,a{t),a{ty, Pit), Pity) = ^ Hi{cu,a{t),a{t)*, P{t), P{t)*) e''' (149) 



and restrict (|146|) - ( |148| ) to these solutions by substituting S = H. This determines the 



autonomous four- dimensional system: 

a = \a - -nK {a* + P+ < g,H^i>) <i),H2> (150) 
(3 = \p -nK{a + P*+ < g,Hi>) < K-'ip,H^2> ■ (151) 

For this bifurcation the mode amplitudes transform according to 

(j)-{a,p) = {e-'^a,e'^p) (152) 
K-{a,p) = {p,a), (153) 

and the 0(2)-invariance of the center manifold implies that the Fourier coefficients of H have 
the form 

Hi = j2 W*y~' hI'\u;, lap, \p\\ aP, a*P*) (154) 

j=0 

for / > 0. Here H^''^ are functions of aP,a*P* and uj that satisfy 

hI^\-uj, lap, |/3|2, aP, a*p*y = Hf~'\u}, |/5p, |ap, a/3, a*/?*) (155) 
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for j = 1,2,...,/. The analysis leading to ( |154| ) - ( |155| ) is analogous to the argument 



given above for steady-state bifurcation: the invariance of the manifold requires the Fourier 
components of H to satisfy 

e-'^'^Hi{uj,a,a*,(3,(3*) = Hi{e,u,e-''^ a,e''^ a* ,e'^ (3,e-'^ (3*), (156) 

and 

Hi{-u, a, «*, /?, /3*y = Hi{uj, A f3*, a, «*), (157) 

The constraint ( |156| ) implies that the function (a*)' Hi is invariant under rotations, hence 
it can only depend on the mode amplitudes only through the basic rotation invariants 
a/?, and a*(3* and must contain an overall factor of (a*)'. By combining |ap 
and a* (3* there are / + 1 distinct ways to generate this factor so («*)' Hi takes the form 



(a 



^^Hi = J2 I"!"' ia*P*y'' hI'\uj, \p\^, a/3, a*p*) (15^ 



j=0 



where the functions ifp'' are arbitrary functions of jap, a(3,a*(3* and uu. Thus Hi has 
the asserted form ( |154| ); the condition imposed by reflection invariance ( |157D implies ( |155| ). 



For / = 1 and 1 = 2, the components needed in ( |150| ) and ( |151| ), we have 



Hi = aH[^\u,\a\^,\p\^,aP,a*p*)+P* H[°\uj,\a\'^,\P\'^,ap,a*P*) (159) 
H2 = a^H^^\uj,\a\^,\(3\\a(3,a*f3*) + a(3*H^^\uj,\a\^,\f]\^,a(3,a*f3*) + 

(/?*)2i7f (u;,|«|M/?|2,a/3,a*/?*). (160) 

Since the center manifold geometry requires that the Taylor expansion of H begins at second 
order, the leading order terms in these components are 

Hi = 0{3) (161) 
H2 = a'/if (cu) + a/5* /i^'^(cu) + (/?*)' /if (cu) + C(3); (162) 
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the third order terms will not be required to determine the cubic terms in the amplitude 
equations. Note that the general relation ( |155| ) implies 

h^^\-u;r = hPiLo) (163) 
hS\-u;r = h^\oo). (164) 



The four-dimensional vector field ( |150|) - (|151| ) commutes with the 0(2) symmetry (|152| ) 



(|153|) on the center manifold; in addition this system can be further simplified by making 



a nonlinear change of variables {a, 13) {a', (3') to put ( |150| )-( p!5TD in Poincare-Birkhoff 
normal form. [p!7| , This change of variables introduces an additional symmetry defined 
by (a,/?) — > e^^{a,(3); this is the "phase shift symmetry" that characterizes the Hopf normal 
form even in bifurcations without symmetry. |pO|, ^ Thus the normal form for ( |151| ) will 
have 0(2) x symmetry and must be of the general form[|l^ |18| 



^1 ^ = [p{u, A) + zq{u, A)] ^ + [r(n, A) + zs{u, A)]6 ^ , (165) 

where 6 = — and p, q, r, and s are real- valued functions of the two basic invariants 
u = + and A = 5^. A final reduction to two dimensions is accomplished by 
introducing amplitudes and phases a' = Xie^'^^ and [3' = X2e^"'^ into ( |165| ) and noting that 
the two-dimensional amplitude equations, 

xi = [p + 5r]xi (166) 
^2 = [v ~ 5r]x2 (167) 

are decoupled from the phase evolution ai = q + 5s and d-2 = q — 5s. This decoupling is a 
consequence of the phase shift symmetry. Upon expanding ( |166| ) and (|167|) we obtain the 
leading order behavior 

xi = [p(0) + (p„(0) - r(0))a;2 + (p„(0) + r(0))x2 + ■ ■ ■]xi (168) 

^2 = [p(0) + (p„(0) + r{0))xj + (p„(0) - r{0))xl + ■ ■ ■]x2; (169) 
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where p{0) denotes the function p{u, A) evakiated at the origin. From the hnear terms in 
( |151| ) we find p(0) = ReA and the third order coefficients p«(0) and r(0) are evaluated in the 
following section. 

Assuming the non-degeneracy condition, r(0)p(0) (Pm(0) — r(0)) 7^ 0, then the amplitude 
equations ( |168|) - (|169|) may be truncated at third order. Jl^ In addition to the incoherent 
state (xi, X2) = (0, 0), setting Xi = X2 = yields two other types of solutions summarized in 
Table 1. There are travelling wave (TW) solutions (xtw, 0) and (0, xtw) with ReX + {pu{0) — 
r(0)) x^w = cind standing wave solutions (SW) (xsw, a;sw) with ReA + 2pu{0) Xg-y^ = 0. 

In the rotating frame, both TW solutions have same frequency cjtw = Q + x'^w^ but 
the two states propagate oppositely in phase; the (xtw, 0) state corresponds to a density 
depending on {6 + u!T^^t) and the (0,xtw) state evolves according to {9 — u!^^^t). Note that 
q and s are evaluated on the TW solution, and therefore near the bifurcation cutw is fixed 
by the eigenvalue (|9TD at criticality {K = Kh). 

^TW — + higher order terms in x^w 



/uj'o-{e + D)^. (170) 

In the lab frame (^), the frequencies of these two states are split to cJ — u^w and U + cjtW) 
respectively. 

The SW frequency is usw = 1 where q is evaluated on the SW solution; near criticality 
tUg^ is also given by ( |17CI| ). The physical appearance of the SW can be described from its 
form near onset. 



psw{0,uj,t) = ^ + KxswIt^ 
Zn 



LJ + -iD Lit - + iD 

+ higher order terms in xsw i^"^^) 

where ao and Po are initial phases. There are two counter-propagating synchronized clumps 
containing oscillators with native frequencies near —uj^-^ and +uj^'w, respectively. Thus 
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Table 1. Branching equations and eigenvalues for non-degenerate Hopf bifurcation 



Solution type (xi, a;2) Branching Equations Eigenvalues 

Travelling wave (TW) ReA + {pu — r) x^w = '^{Vu — f) a^xw 
(xtw, 0) and (0, xtw) ReA + {pu + r) x\-^ 



Standing wave (SW) ReA + 2p„ Xg-^y = 4p„ x 

(a;sw,a;sw) -Ar x\ 



2 

SW' 
2 

SW 



NOTES: 

The coefficients Pu and r are evaluated at (xi,X2) = (0,0). 

the population exhibits two simultaneous macroscopic oscillations. In the lab frame these 
macroscopic oscillators appear as a slow oscillator at frequency uJ — wsw and a fast oscillator 
at cU + wsw- 

The direction of bifurcation for these waves is determined by (Pm(0) — r(0)) and Pm(0). 
Their linear stabilities are easily calculated: the TW eigenvalues are 2(pu(0) — r(0)) and 
ReA + (Pm(0) +r(0))a;^w; the SW eigenvalues are 4p„(0)a;g^ and — 4r(0)a;g^. Depending 
on the signs and relative magnitudes of Pu(0) and r(0) there are six possible bifurcations as 
shown in Fig. 2. For the bimodal population ( |53|) we shall see that (pn(0) — r(0)) and Pu(0) 
are negative while r(0) is positive. Thus the SW are supercritical and stable and the TW 
are supercritical and unstable. 

3.2.2 Evaluation of p^ and r 



From ([T62D for H2, 

<tp,H2> = < ^, h^^^ > +af3* < h^^^ > +{(3*f < tP, hf > +0(3) (172) 
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(3^ <K-i,, h^^^* > +C(3) 



(173) 



which determines the nonhnear terms in ( |150| )-( p!51| ) through third order: 



a 



a 



a[< ij, > < i,, h'i^ > 



7tK 



(1) 



P[< K ■ ^, /if* > <K-ij, h^2^* > \a\ 



a*[af3* < ^, h^2^ > < ^, /if >] 



-7fK 



(0) 



(1)* 



\ a[{a*f < ^ ■ ^, >+P^<^.^p, hf* >] ) + ^^^^ ^^^^^ 



Only the first set of nonhnear terms in (|174|) have the phase shift symmetry and are 
retained in the normal form: 



a' 



A 



a 



the remaining cubic terms and all fourth order terms in ( |174|) are removed by the normal 
form change of variables (a,/?) — ^ {a', (3'). From ( |163| ) and ( |164| ) the coefficients in ( |175| ) 
are related: < ijj, h2'^ >=< k ■ -ijj, /if* > and < -ijj, /if >=< k ■ 'i/j, /if* > so that the result 
in ( |175| ) can be rewritten in the standard form ( |165| ). Noting that = {u — 6)/2 and 
IP'l"^ = {u + 5)/2, ([T7|) becomes 

' nK 



a' 



X- 



(< ip, /if > + < V, h''2' >)u + C(4) 



(1) 



+ 



ttK 



K^l^.hf >-<^,h'i> >) + 0{2) 



(1) 



a 

J' 
a' 



(176) 
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By comparing ( |165| ) and ( |176| ) we can identify the cubic coefficients pu and 



Pu{0) 

r(0) 



nK 
nK 



Re < ^, [/i^^^ + h^2'] > 



,(2)l 



Re < tfj, [h^2^ - h'^^'] > 



,(2)i 



(177) 
(178) 



in terms of the center manifold coefficients h^^?^ and h^o\ 



For this bifurcation the equation for H ( |127| ) becomes 



dH . OH ■ 



(179) 



CH - TcK [e'\a* + P+ < g,H^i >)[H2- <ij,H2> ^- <K-i)\H2 > k ■ i)*] 
+2e'^'[{a*+p+<g,H_i >)Hs{u) ~ {a + P*+ < g, >){aij + p* k ■ ^j* + Hi)] 
+ ^ / e''' [{a* + f3+<g, H^i >)Hi+i{u) - {a + (3*+ < g, >)Hi_i{uj)] + cc} , 



l>3 



and the functions /i2^\ h2^\ and /i2°^ may be determined from (|179|) . By inspection of the 
/ = 2 component we note that the inhomogeneous terms 2tt K e^"^^ [a + (3*) {aip + (3* k ■ tp*) 
in ( |179| ) generate the leading terms /i2^(c<j) in H2; since Hi{uj) ~ 0{3) for / 7^ 2, applying 
( |162|) to the / = 2 component of ( |1 79|) and equating terms of the same order, we obtain 



2A h^2^ 
(A + A*) h^^'> 
2X* 



L2h^2^ + 27iKi; 

L2h^2^ + 2n K + K ■ ^*] 

h2h^2^ + 27: K K ■ ^\ 



;i80) 
;i8i) 
;i82) 



with solutions 



—in K ipioj) 
uj — Zq — i2D 
-i7f K (ipiuj) + K ■ ipiu)*) 

UJ — i Im Zq — i2D 
-in K K, ■ 

UJ + z*r,- i2D 



(183) 
(184) 
(185) 
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These results satisfy the relations ( |163| ) - ( |164| ). 

From ( |177] ) - ( |178| ) these results determine the cubic coefficients explicitly 



r(0) 



Im 



Im 



+ 



{uj — ihxv zq — i2D) ' {uj — zq — i2D) 



;i86) 



;i87) 



[uo — ilm zq — i2D) {oo — zq — i2D) 

Provided Im zq > —D, any singularities in the lower half plane are due to g{uj). Hence 
by inspection we can see that at criticality for the bifurcation (Im zq = 0), the limit D 
is finite for both Pu{0) and r(0). By contrast the corresponding limit for the Vlasov 
equation involves pinching singularities that produce a singular behavior. [ 16 1 Here the 
well-behaved limit a.t D = suggests that the oscillations produced by the bifurcation will 



grow like J\K — Kc\ even in the absence of noise. 



3.2.3 Example of a bimodal population 



We illustrate these results for the bimodal population described in (BH); for Auo > K the 
eigenvalues (A, A*) are given by A 



-izQ = {K-Kh)/4:±inwheTe Q = - {K/Af and 
Kh = 4(e + D). UK increases above Kh, there is a Hopf bifurcation. The coefficients Pu{0) 
and r(0) can be evaluated analytically at the bifurcation K = Kh, by closing the integration 
contours in the lower half plane and using fi^ = cJq — (e + D)"^ one obtains 

iKh + {n - i{e + D)y 



2 [oo^^ - {Q - t{e + D)y] 



(188) 



Re < ip, h^2^ > 



Re < ip, hf > 



IQixKhD (e + Dfn'^ujl 



|a;o + i(e + 2D)|2 \uj^^ - {^ + ^{e + D)Y P 



A'KKhn'' (e + 2D) ul [Aul + (e + 2D)f - (e + D)f 
\ujI + (fi - i(e + L'))2|2 \ujI -{n- i{e + 2D)y\^ 



(189) 



(190) 
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The direction of bifurcation and stability of the travelhng waves and standing waves are 
controlled by Pu{0), r{0), and Pu{0) — r{0). Since Re < ijj,h^2^ > and Re < ijj,h^2^ > are 
clearly positive, one can see by inspection that both p„(0) and Pu{0) — r(0) are negative, and 
the SW and TW solutions are always supercritical. The bifurcation produces stable TW 
solutions if r(0) < and stable SW if r(0) > 0; see Fig. 2. Subtracting Re < ip,h2 > from 
Re < V^, h2^^ > gives 

\ujI - {Q + i{e + D)y \^ lui - {Q - i{e + 2D)y\^ looo + ^{e + 2D)\^ ^ ' 

where 

P{e,D,uj^) = 4(e + 2D)cu^ + (e + 2D)[4e^ + 14eL> + llDVo 

+D{e + D) (3e + 3D) [e^ + htD + (192) 

Thus for this bimodal population r(0) is always positive and the bifurcation always leads to 
stable standing waves. 

In light of Kuramoto's earlier discussion of the partially synchronized bimodal population 
(cf. the remarks in Section 2.3), it is interesting to compare the frequencies Wpeak in (§§) and 
'^sw (|170|) for the SW near onset in the relevant frequency range (l + Z^/e) <co'o/e<oo. 
From dH) and (|l7g) we find 

V'^P«.t/ [(2w„/e) - \/l + (wo/e)"] ^1 + ('^o/e)' 

as the peaks move far apart {uj^/e —>■ oo) this ratio approaches unity and the oscillators in 
the peaks are indeed the oscillators that synchronize. However as c^o/e approaches (l + D/e) 
the oscillators comprising the standing wave synchronize at frequencies that are substantially 
shifted away from the peaks of the native distribution towards the interior of the distribution; 
this is illustrated in Fig. 3. 
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3.2.4 Comparison to Bonilla, Neu, and Spigler 

Bonilla, et al. have previously calculated the amplitude equation for a Hopf bifurcation 
from the incoherent state, but they neglected to include the second eigenvector k ■ \1/ in ( ^7|) 
that is forced by 0(2) symmetry. In our notation this omission has the effect of setting the 
amplitude (3{t) of this mode to zero; if this is done then the two-dimensional system for 
(xi,X2) in (|168|) - ( p.69| ) is reduced to0 

xi = [p(0) + (p„(0) - r{0))xl + ■ ■ ■]xi. (194) 

Since Pu(0) — r(0) < for the bimodal population, this equation predicts a single TW 
(xi, X2) = (xtw, 0) which is supercritical and stable. The SW solutions (xi, X2) = {xsw, a^sw) 
and the second TW solution (xi, X2) = (0, xtw) are omitted; in addition the actual instability 
of the TW is missed because the TW eigenvalue ReA + (pu(0) + r(0)) x^^^ corresponding to 
perturbations in the direction of the SW state is not included. 

3.2.5 Comparison to Okuda and Kuramoto 

Okuda and Kuramoto consider two populations of oscillators {i = 1,2) with number 
densities nj(0, t); in the absence of couplings each population has a definite native frequency: 
uji = uj -\- Auj/2 and UJ2 = uJ — Auj/2 respectively. The evolution of ni{(j),t) and n2{4>,t) is 
described by 

"'^i^T + 'al 1 ^'i^i^' / sin(0 - 0') [K[ni{(f)', t) + K'n2{(p', t)] 



dt dcj) (90 I Jo 

= ~^2-^ + 1^2(0, t) d(j)' sin(0 - 0') [^2^2(0', t) + K'ni{(t)\ t)] 

+D^. (196) 



'^This corresponds to equation (3.18a) in Bonilla, et al. with the identification Kxi/2 — > \A\. 
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The intra-population couplings are K[ and K2 and the couphng between the populations 
is K' . The incoherent state is an equilibrium described by ^i(0) = "^2(0) = l/^^r, and the 
authors analyzed bifurcations from this state. This model has 0(2) symmetry if a; = and 
K[ = K^. 

For the bimodal population (^) at e = 0, the Fokker-Planck equation (^) reduces to 
the 0-K model but with the constraint of equal couplings K[ = K2 = K' . For the e = 
distribution, 

(7(0;) = i[5(cu + L^o) + 5(^-cuo)], (197) 

we define number densities by integrating over the frequency distribution of each sub- 
population 



ni{e,t) = I duj6{uj - LUo) p{e,uj,t) (198) 

(199) 

n2{e,t) = I du6{uj + uJo) p{e,uj,t) (200) 



so that M) becomes 



dni drii d 



-UJo- 



dt " 89 

dn2 dn-i d ( f'^^ 



+ ^ [ni{e, t) J^^ dO' sin(^ - e') [{KI2) n,{e\ t) + (ir/2) n^iO' , t)] 



(201) 



+Uo- 



dt 06 



+ ^ [n2ie, t) " de' sm{e - e') [{K/2) n^ie', t) + {K/2) n,{e\ t)] 



(202) 



The equations (|9|) - (|T9|) have this form if we set Aw = 2a;o and K[ = K'^ = K' = K/2 
and transform to the rotating frame 9 = cf) — uJt where average frequency is zero. One of the 
cases considered by Okuda and Kuramoto is K[ = K'2 = 1 with K' and Au variable. Their 
results for the steady state bifurcation and Hopf bifurcation are consistent with those of this 
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paper. In particular, for the Hopf bifurcation they find that the SW solutions are expected 
to be stable. 

3.2.6 Perturbing 0(2) ^ S0(2) 

If the distribution g{uj) is perturbed so that the reflection symmetry (|l^) is broken, then 
the complex eigenvalues of multiplicity-two split into two non-degenerate complex conjugate 
pairs. The Hopf bifurcation considered here perturbs to two distinct Hopf bifurcations cor- 
responding to the distinct pairs of eigenvalues. The resulting Hopf bifurcations individually 
lead to time-periodic states in the form of rotating or travelling waves and the interac- 



tion between these waves yields modulated waves via secondary bifurcations, [p^j - p7| The 



modulated waves are the perturbed form of the standing waves discussed above; Knobloch 



discusses such perturbed bifurcations in more detail. fT^ 

4 Discussion 

The bifurcation analysis is summarized in Fig. 4 which shows the stable nonlinear states 
that arise from instabilities of the incoherent state. Our discussion does not treat the transi- 
tion from the stable partially synchronized state to the stable standing waves; this transition 
is represented by a dashed boundary in Fig. 4. The interaction of the standing waves and 
the synchronized state can be studied by analyzing the co dimension- two bifurcation defined 
by the merger of the Hopf and steady state bifurcations. The linear theory of the non- 
semisimple real eigenvalue was developed in Section 2 of this paper and a normal form theory 
for the 0(2) Takens-Bogdanov bifurcation has been given by Dangelmayr and Knobloch [p8 



One complicating feature of the bimodal population is the subcriticality of the steady state 
bifurcation at the codimension-two point; one would like to include the stable partially syn- 
chronized state at finite amplitude in the analysis. This could be done by generalizing the 
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Dangelmayr-Knobloch theory to combine the nonhnear degeneracy of the steady state bi- 
furcation with the hnear degeneracy. Formally this leads to a codimension-three bifurcation 
and would require a normal form of higher order; the analysis of Dangelmayr and Knobloch 
truncates the normal form equations after cubic terms. 
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6 Figure captions 

1. Eigenvalues of the incoherent state as a function of e,D, K, and uoq. The solid curves 
indicate parameter values at which there is an eigenvalue on the imaginary axis; along 
K — Kc there is an eigenvalue at zero and along K — Kh there is a pure imaginary 
complex-conjugate pair. At the juncture K = 4(e + D) and utq = {e + D), is the 
codimension two point where the Hopf bifurcation surface intersects the steady state 
bifurcation surface. The dashed curves indicate parameter values where eigenvalues 
coalesce with the continuum; upon crossing K = K'^ or K = Kg an eigenvalue either 
appears or disappears. Along the dot-dashed line K — Aluq, the two eigenvalues collide 
on the real axis. The inset diagrams show the qualitative features of the eigenvalue 
spectrum in each region. The continuous spectrum at —D is not drawn. 

2. Bifurcation diagrams for Hopf bifurcation with 0(2) symmetry. The stability and 
direction of bifurcation for the travelling waves (TW) and standing waves (SW) are 
determined by the two cubic coefficients r(0) and Pu(0). The incoherent state is the 
horizontal branch and the two TW solutions are drawn as a single branch. Stable 
solutions are solid branches and unstable solutions are dashed. The signs of the real 
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parts of the eigenvalues are indicated; see Table 1 for the TW and SW eigenvalues. 

3. Comparison of the standing wave frequency cu^^ and o^peak for the bimodal distribution 
at criticality K = Kh for the Hopf bifurcation along the line {1 + D/e) < (o^o/e) < oo. 
Also shown is the ratio of ujq to ujpeak', the curves are drawn for D/e = 1. 

4. Nonlinear states for the bimodal distribution: the incoherent state (I), the partially 
synchronized state (PS), and the standing waves (SW). The solid lines are the stability 
boundaries for steady state bifurcation and Hopf bifurcation. The boundary between 
the standing waves and the partially synchronized state is shown schematically as a 
dashed line; the precise nature and location of this boundary has not been determined. 
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